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Abstract. We present Monte Carlo simulations of a coarse-grained model 
for Langmuir monolayers of amphiphile molecules on a polar substrate. The 
molecules are modelled as chains of Lennard- Jones beads, with one slightly larger 
end bead confined in a planar surface. They are simulated in continuous space 
under conditions of constant pressure, using a simulation box of variable size 
and shape. The model exhibits a disordered phase (corresponding to the liquid 
expanded phase), and various ordered phases (corresponding to the condensed 
phases) with different types of tilt. We calculate the phase diagrams and charac- 
terize the different phases and phase transitions. The effect of varying the chain 
stiffness is also discussed. 



I. INTRODUCTION 

Monolayers of amphiphiles at surfaces (Langmuir monolayers) have attracted longstand- 
ing scientific interest for various reasons fl|-|3[]: Surface properties of materials can be mod- 
ified and taylored by coating the surfaces with amphiphiles. Langmuir monolayers can be 
exploited to engineer thin film materials with well-defined structures on a molecular level. 
On the other hand, lipid monolayers on water are experimentally fairly accessible model 
systems for biological membranes. Last not least, Langmuir monolayers are experimental 
realizations of two-dimensional systems, which allow to study ordering phenomena in low 
dimensions. 

Experimentally, Langmuir monolayers have been investigated for a long time by mea- 
surements of pressure-area isotherms M. More recently, a number of powerful microscopy 
techniques have been developed, such as fluorescence microscopy and Brewster angle mi- 
croscopy, which have provided insight into the mesoscopic structures in monolayers. The 
emerging pictures for monolayers on water is qualitatively similar for phospholipids, long 
chain alcohols and esters: At low surface coverage, the molecules hardly interact with each 
other and build the two-dimensional equivalent of a "gas". Upon compression, a first or- 
der transition to a fluid-like "liquid expanded" (LE) phase is encountered, followed at even 
higher surface coverage by a second discontinuous transition into a "liquid condensed" (LC) 
area. The transition from liquid expanded to liquid condensed has an important equivalent 
in bilayers, the "main transition", which may be biologically relevant, since it takes place 
at temperatures close to the body temperature for some of the common phospholipids. The 
condensed region contains a variety of different phases, characterized by different types of 
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ordering, i.e., collective tilt order of the hydrocarbon chains, orientational order of the back- 
bones of the chains, and crystalline positional order. A generic phase diagram for fatty acid 
monolayers is shown in Fig. [I] ||. The lowest density phases which coexist with the LE phase 
are typically hexatic rotator phases, i.e., the backbones rotate freely around, and positional 
correlations decay exponentially, but the directions of nearest neighbors are nevertheless 
well-defined. 

Theoretical treatments of Langmuir monolayers have followed three different lines. On 
the one hand, phenomenological descriptions of the different condensed phases in terms of 
Landau expansions in the characteristic order parameters [HH have offered valuable insight 
into the nature and the interrelations of different phase transitions on a very general level. 
On the other hand, Molecular dynamics simulations of atomically realistic models have com- 
plemented experiments and provided structural information on quantities, which are hard 
to access experimentally These two approaches are in a sense antipodal: Whereas 

phenomenological treatments focus on universal properties and make little or no contact to 
the microscopic structure of the systems, atomically realistic models seek to imitate nature 
as faithfully as possible, and to reach quantitative agreement. Hence they account for many 
more details than are actually needed to produce a certain phase behavior, rely heavily on 
the availability of good force fields, and their study is computationally costly. 

As a third line of approach, idealized microscopic models are constructed which incor- 
porate only a few properties of a material, believed to be essential for a given behavior. 
Thus they bridge between phenomenological and realistic models, and relate microscopic 
and macroscopic quantities in a qualitative and semi- quantitative way. 

The question, which features of amphiphiles are essential in Langmuir monolayers, can 
of course not be answered universally. It depends on the region in phase space one wishes 
to study. Attractive interactions between the amphiphiles are important for most phase 
transitions. As long as one studies condensed phases, it is often sufficient to model the 



amphiphiles as anisotropic stiff objects. Grafted rigid rods exhibit tilt transitions fl2Hl4 



molecules with non-circular cross-sections show rotator transitions ||15|| . For the transition 
between the liquid condensed and the liquid expanded phase, however, the conformational 
degrees of freedom of the chains play a crucial role |T0|JT6| , p!7[| . They have been incorporated 
in a heuristic way as "internal degeneracies" in Ising-type two-dimensional lattice models 
for monolayers and bilayers, e.g., in the Pink model fll8,19fl. The interdependence of chain 
conformations and effective chain interactions has to be put in by hand in this approach, and 
a large number of input parameters is required. Models which aim to study more directly 
the interplay of chain conformations and phase behavior have to retain the chain character 
of the amphiphiles explicitly. 

A suitable idealized model for Langmuir monolayers thus represents the amphiphiles by 
flexible chains of mutually attracting monomers, which are grafted to a surface at one end 
("head"). Such models have been formulated on the lattice [20-EM| and in continuous space 

EMI- 

Lattice models can be simulated more efficiently than off-lattice models, yet they can 
produce rather awkward lattice effects especially when orientational order (tilt order) comes 
into play ||24j| . An off-lattice bead-spring model of Lennard- Jones beads has been studied by 



Haas et al p5|[26f and by us p7| under constant volume and constant pressure conditions. It 
was found to display a tilted and an untilted phase, in which the chains are basically arranged 
on a (possibly distorted) hexagonal lattice, and a "fluidized" phase which is reminiscent of 



2 



the liquid expanded phase. Hence it seems a promising candidate for a minimal model, which 
contains only the basic elements responsible for the main transition in Langmuir monolayers. 
Nevertheless, no systematic study of the phase behavior has been presented so far. 

This is the objective of the present paper. We have performed Monte Carlo simulations 
of a bead-spring model very similar to the one used by Haas et al. The models only differ 
in the treatment of the heads: Whereas the head beads in Haas et al's version are identical 
with the chain beads, our heads are slightly larger. We chose this variant in order to ensure 
that the dominant reason for chain tilting in our model is similar to the most common one 
in nature: Tilt is induced by the mismatch of head and tail size. In the model of Haas et al, 
the chains tilt, because they can then "hook" into each other and thus pack more efficiently 
The details of the tilt order (tilt angle, tilt direction etc.) result from a complicated interplay 



between monomer packing and chain stretching ||26|| , which is highly model dependent and 
has probably little to do with the factors which influence the tilt in real monolayers. On 
simple geometrical grounds, two of us have argued earlier that the direction of tilt depends 



on the size of the head groups |30| . There is also experimental evidence for such a connection 



With our choice of the head size, we ensure that the model exhibits two different tilted 
phases at zero temperature, a low-pressure one with tilt towards nearest neighbors, and a 
higher-pressure one with tilt towards next nearest neighbors. 

Our paper is organized as follows: In the next section, we specify the model and comment 
on some aspects of the simulation techniques and the data analysis. The results are presented 
in section 3: We characterize the phases and phase transitions, show the phase diagrams, 
and discuss the effect of the chain stiffness. We summarize and conclude in section 4. 



II. MODEL AND TECHNICAL DETAILS 



Following Haas et al [25,2q|, we model the amphiphiles as chains of beads, which are 



connected by springs of length d subject to the spring potential 

, -!fd 2 s \n(l - (d-d ) 2 /ds 2 ) for \d - d \ < d s 
V s (d) = { 2 V ^ (1) 

oo for \d — d \ > ds 

This so-called "Finite extension nonlinear elastic" potential (FENE) is basically harmonic 
at d ~ do and has a logarithmic cutoff at d = do ± d$. Furthermore, we impose a stiffness 
potential 

V A = k A - (l-cos0) (2) 

on the angle 9 between subsequent springs. The stiffness potential favors angles 8 = 0, i.e., 
straight chains. Beads are not allowed to enter the half space z < 0; moreover, one end bead 
of each chain (the "head") is confined to remain within the plane z = 0. Thus we assume 
a very strong binding force between the hydrophilic head group and the water surface, and 
the latter is approximated by a perfectly sharp and flat interface. Tail beads interact via a 
truncated Lennard- Jones potential 
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, . air) 12 — 2(alrf + v c ) for r < 2a 
V L j{r) = { vv 7 > ~ , (3) 

for r > 2a 

where v c = 127/4096 ~ 0.031 is chosen such that Vlj(t) is continuous at r = 2a. The 
interactions between head beads are purely repulsive, 

\e H - ((a H /r) 12 - 2(a H /rf + l) for r < a H 
V H (r) = < v 7 . (4) 

for r > a H 

The attractive part here has been cut off for reasons of computational efficiency. Note that 
the head size an differs from the tail bead size a. Head and tail beads interact with a 
repulsive potential of the form (||), in which an is replaced by (an + c)/2. 

The parameters e and a define the units of energy and length. To complete the definition 
of the model, we have to specify the remaining parameters do, ds, ks, kA, £h, and an- Our 
choice was motivated by the idea that one bead should represent roughly two CH 2 groups 
in an actual alkane chain. Comparing a straight model chain with an ideal all-trans state 
hydrocarbon chain, with realistic potential parameters of united-atom potentials taken from 
the literature (e.g., from Ref. [B3]), one finds that the bond length d Q should be approximately 
0.7 times the chain diameter, do = 0.7a. The identification also allows for a rough estimate 
of the absolute values of a and e: a ~ 3.8A and e ~ 240ksK, where ks is the Boltzmann 
constant. These values should of course not be taken too literally, since the model is much 
too simple to allow for quantitative comparisons with experimental systems. 

The spring constant ks was chosen very strong, ks = lOOe, such that the lengths of the 
springs are approximately constant at all temperatures of interest. The value of the cutoff 
ds then has little influence on the properties of the system; we use ds = 0.2a. The stiffness 
constant kA can be estimated by adjusting the average (cos 9) of a single free chain in our 
model at a given temperature to the corresponding value in a single free alkane chain. Such 



an estimate would yield kA ~ 5e at room temperature. Haas et al |25| , |26|| have used kA = 10e. 
Here, we have mostly used the same value (kA = 10e) in order to be consistent with their 
work. For the reasons mentioned in the introduction, the size of the head beads was taken 
to be ajj = 1.1a. The influence of the head size on the phase behavior shall be discussed in 



detail elsewhere [p9| , p3| . The prefactor e# was chosen = e. 

The simulations were performed at constant spreading pressure in a simulation box of 
variable size and shape. More specifically, we study n chains of length N on a parallelogram 
with side length L x and L y and angle a. Periodic boundary conditions were applied in these 
two directions, and free boundary conditions in the third. Our Monte Carlo moves include: 

• Attempts to displace single beads 

• Attempts to vary L x , L y or a, i.e., to rescale all coordinates such the the configuration 
is stretched or squeezed in one direction, or sheared ( "volume moves" ) 

The trial moves are accepted or rejected according to a standard Metropolis prescription 
with the effective Hamiltonian 



H = E + UA — nNT ln(A), (5) 
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where E is the internal energy, II the applied spreading pressure, and A = L x L y sin a the 
area of the simulation box. We have also implemented collective moves, in which chains were 
displaced as a whole, and volume moves, in which only the coordinates of the head beads 
are rescaled, but inner molecular distances and angles are kept constant. The Hamiltonian 
(0) then has to be replaced by 

H = E + UA-nT\n(A). (6) 

Unfortunately, these collective moves did not reduce the time needed to generate uncorre- 
lated configurations significantly. Similarly, we have implemented continuous configurational 
biased Monte Carlo moves ||35|| , but found that they brought no improvement in our partic- 



ular system. 

In order to check that no internal stress is present in our simulations, we have determined 
the internal pressure tensor 



where the sum i runs over all monomers, a, (3 over the x and y coordinate, denotes the 
force acting on monomer i, and 5 a p is the unit matrix. According to the virial theorem, 
11^ should be diagonal and identical to Ii5 a p at mechanical equilibrium. This was the case 
in our simulations, if we used a simulation box of variable shape. In simulation runs with 
a rectangular box, we sometimes obtained nonzero off-diagonal elements U xy in the tilted 
phases. 

We will present results for n = 144 chains of length N = 7. The average decorrelation 
time lies between 200 and 1000 Monte Carlo steps (MCS), where one MCS consists of 
Nn = 1008 attempts of monomer moves, and one attempt to rescale L x , L y , and a. In 
general, the systems were equilibrated during 70.000 MCS, and data were then collected 
from every 500st configuration over a period of at least 200.000 MCS. 

The simulations were supplemented by a low temperature analysis. The zero temperature 
ground state was determined by minimization of the enthalpy (|5]). A harmonic expansion 
was then performed in order to determine the free energy G at some given low (nonzero) 
temperature To. Given this reference value, one can calculate the free energy at other 
temperatures and pressures from simulations by means of a thermodynamic integration 

g(ilt) = G(n ,T ) + k B r j r {dn'-^j - dT 'j^}^ (») 

as long as the path T from (n ,T ) to (II, T) does not cross a first order phase transition. 
By comparing the free energies of different states, we have localized the transition points 
between phases at low temperatures where hysteresis effects were strong. 



III. RESULTS 

Figure § shows temperature-area isobars for a selection of low pressures (a) and high 
pressures (b). One clearly observes a jump in the area per molecule, which moves to higher 
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temperatures as the pressure increases. At high pressures, one discernes in addition a kink 
at low temperatures, indicating the presence of a second phase transition. 

The phases can be characterized by the typical features of the pair correlation functions 
(Figs. ||and §) and structure functions (Figs. | and||). For example, the two dimensional 
correlation functions in the intermediate temperature state at high pressures are precisely 
those of a hexagonal lattice. Fig. |3] shows pair correlation functions for the head groups, 
the projection of center of gravity of the chains onto the xy plane, and the points where the 
chains pass through the plane at z = 2a above the surface at pressure II = 100e/a 2 and 
temperature T = le/ks, which is slightly above the first phase transition. The three curves 
do not differ from each other qualitatively, and the position and relative heights of the peaks 
are consistent with those of a hexagonal structure. At temperatures below the first phase 
transition or at lower pressures, each of the peaks splits up in two. This indicates that 
the hexagonal lattice is distorted in one of the high symmetry directions, either the nearest 
neighbor or the next nearest neighbor direction (for intermediate directions the peaks would 
split up in three). An example is shown in Fig. || (see the curves for the lowest temperature 
T = OAe/ke)- From the large height difference of the twin peaks, one can infer that the 
lattice is stretched in the direction of nearest neighbors in this specific case. The direct 
inspection of configuration snapshots reveals, not surprisingly, that the lattice distortion 
goes along with a collective tilt of the chains in the direction of the distortion. At low 
pressures (II ~ 10e/cr 2 ), the hexagonal lattice in the tilted phases is stretched by roughly 
10%. 

With increasing temperature, the structure of the correlation functions is gradually lost. 
Slightly below the phase transition, the correlation functions of the head lattice are fluid-like, 
with peaks of monotonically decreasing height for the first, second and third coordination 
shell. They do not change qualitatively as the phase transition is crossed (Fig. £| (a)). 
In contrast, the correlation function for the projections of the center of gravity still shows 
some solid-like structure right below the phase transition, and loses almost every structure 
right above the phase transition (Fig. ^ (b)). In the high temperature state, the head 
positions are much more correlated than the chain positions. We conclude that the phase 
transition associated to the area jump is a melting transition, and that it is driven by 
the chains. The chains maintain the order below the transition, and promote the disorder 
above the transition. This is consistent with results from molecular dynamics simulation by 



Karaborni and Toxvaerd HTT| of a realistic model. 



The structure function is defined by 



^ nN 



^exp(igf.,) , (9) 



where the sum runs over all monomers in the system. Note that in a finite simulation 
box with periodic boundary conditions, S(q) for a specific configuration is only defined for 
vectors q whose projection on the xy plane are sums of integer multiples of the basis vectors 

and 

However, the dimensions of the box fluctuate in our simulations, hence the basis vectors 
fluctuate as well. In order to overcome this problem, we have laid a fine-meshed grid on the 
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xy plane and summed up all the contributions to S(q) within a mesh. Fig. [| (a) and (b) 
shows the resulting structure factors in the plane of q z = for an disordered state (a) and 
an untilted ordered state (b). The structure factor of the disordered state is isotropic and 
shows the usual features of a fluid structure factor. In the untilted ordered state, one finds 
the Bragg rods of the hexagonal lattice. They are sharply peaked in the xy plane, but have 
a considerable width in the z direction, thus the term "rods". In the tilted ordered state, 
the plane of maxima tilts such that it stays perpendicular to the long axis of the chains [[3] . 
Thus the peaks belonging to q vectors which are not perpendicular to the tilt direction move 
out of the q z = plane. This is illustrated in Fig |]for a state with tilt towards next nearest 
neighbors. The internal structure of the rods in the z-direction reflects the structure of the 
monolayer. For example, the width of the rods is inversely proportional to the width of the 
layer, and every rod is surrounded by a multitude of weak "satellite maxima" which are 
caused by the sharp steps in the density profile at z = and at the outer surface. After six 
low satellite maxima, another strong peak is found, reaching a height comparable to that of 
the main peak. These peaks reflect the "periodic" arrangement of monomers within a chain. 
They are found at distances of approximately Aq z rs 2ft /do cos 8 and integer multiples from 
the main peak, where do is the favored distance between monomers (see eqn. ([[])). Their 
appearance is a very specific property of our simulation model, and not interesting from a 
general point of view. Hence they shall not be studied any further. 

In order to quantify our findings, we have analyzed a number of suitable order parameters. 
For example, we determine the hexagonal order parameter of two dimensional melting 



Here the first sum j runs over all heads of the systems, the second k over the six nearest 
neighbors of j, and <pjk is the angle between the vector connecting the two heads and an 
arbitrary reference axis. The quantity thus measures the orientational long range order 
of nearest neighbor directions. It is nonzero in the hexagonal (quasi) crystalline phase and 
in the hexatic phase. As an order parameter which describes the collective tilt of molecules, 
we have computed 



which corresponds to the length of the average projection of the head-to-end vector of the 
chains on the xy plane. Here [x] and [y] denote the x and y component of the head-to-end 
vector, averaged over the chains of a configuration, and (■) denotes the thermal average 
over all configurations. The quantitity R xy is nonzero in phases which break the azimuthal 
symmetry, i.e.. phases with collective tilt, and zero otherwise. Note that the average tilt 
angle 9 between the head-to-end vector of the chains and the surface normal is always 
nonzero. 

The quantities and R xy are shown as a function of temperature for various pressures 
in Figs. [7] and || The area jump in the isobars goes along with a drop to almost zero of the 
melting order parameter ^>§. This substantiates our earlier speculation that the transition 
corresponds to a melting transition. Furthermore, we infer from the decrease of R xy with 
the temperature, that there is also a tilting transition from a collectively tilted phase at 
small temperatures, to an untilted phase at high temperature. The melting transition and 




3=1 k=l 



(10) 




(11) 
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the tilting transition occur simultaneously at low pressures, and decouple from each other 
at high pressures. The tilting transition then precedes the melting transition and seems to 
be continuous. 

At our small system size, it is not possible to decide whether the ordered phase is crys- 
talline or hexatic. Moreover, we are not able to establish unambiguously the order of the 
melting transition. These are two closely related issues of high interest. Even in much sim- 
pler two dimensional systems (hard disks, Lennard- Jones disks), the question whether they 
melt discontinuously in one stage, or continuously via a hexatic phase in two stages, is 
still a matter of debate. The transition from a hexatic to a fluid phase is usually believed 
to be continuous. In the case of amphiphile monolayers, however, we have argued that it 
can be driven first order, as an effect of the interplay between chain entropy and chain 
packing ||17|| . We have already noted that the melting transition in our system is mainly 



driven by the chains, which enhances the likelihood of such a scenario. The transition may 
also be discontinuous at low and intermediate pressures, and continuous at high pressures. 
The pronounced jumps observed in our simulations seem to indicate a line of discontinuous 
transitions; on the other hand, we have not encountered significant hysteresis effects except 
at very low pressure, II = 1. Simulations of much larger systems and a thorough finite-size 
analysis would presumably be necessary to distinguish between first order and continuous 
transitions. 

It is instructive to also consider the distribution of tilt angles 9. Let us first look at the 
average (9) (Fig. §). In the low pressure regime, where the melting and the tilting phase 
transition coincide, it drops down at the transition and then rises slowly with temperature. 
At higher pressures, where the two transitions decouple, it first decreases with temperature 
until the tilting transition is passed, then stays low in the temperature region of the untilted 
ordered phase, but jumps to a higher value at the melting transition. The jump is related to 
the jump in the area per molecule at that transition: the molecules have more space to lie 
down. The average tilt angle is coupled to the molecular area A/n by the requirement that 
the bead density in the monolayer should not vary much, i.e., the total volume occupied by 
the monolayer is close to constant. In the condensed region, where the chains are mostly 
straight and aligned, this implies that the quantity A cos(6 l ) jn is approximately constant and 
equal to a c , the area per molecule in the untilted high pressure phase. Such a dependence 
has indeed been reported experimentally [j37 l. Similarly, we find here that the product of A 
and cos((6 1 )) depends much less on the temperature and pressure than the area per molecule 
A/n itself (Figure [10]). In particular, its value right below the melting transition is found 
to be a c ~ 0.985cr 2 at all pressures except for the very highest, II = 100e/cx 2 , regardless of 
whether the condensed phase is tilted or untilted. Hence the volume density in the monolayer 
seems to trigger the melting transition rather than the area density - which corroborates 
our earlier assertion that the melting transition is driven by the chains. 

Figure O shows the histogram of the tilt angle P{9)/ sin# at pressure II = 50 for dif- 
ferent temperatures. Below the tilting transition, P{9)/ sin 9 has a clear maximum. As the 
temperature is increased, the maximum moves down towards lower values 9. At the tilting 
transition, it merges into 9 = 0. From there on, it becomes broader, which explains the 
increase of (9) at higher temperatures. 

Finally, we turn to the discussion of the direction of the tilt. It can be determined from 
a histogram of the angle between the momentary tilt direction and the bonds connecting 
nearest neighbors. If the tilt angle is well-defined, this histogram should have six peaks, and 
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their positions indicate the direction of tilt. At low temperatures T ~ 0.1, we find two phases 
with well-defined tilt directions towards nearest neighbors and next nearest neighbors. The 
transition between them is strongly first order, and the thermodynamic integration methods 
described in the previous section had to be used to locate the transition points. At higher 
temperatures, the transition washes out, and in some regions of phase space it is hard to 
determine whether the tilt direction is at all locked to the underlying hexagonal head lattice. 
In order to quantify the "locking", we define an order parameter $6, which is very similar 
to the hexagonal order parameter (eqn. (|10D ). 



j n d 



(12) 



The notation corresponds to that in eqn. (|T0|) , except that is now the angle to the average 
tilt direction in the current configuration rather than simply that to an arbitrary reference 
axis. The crucial difference to the definition of ipe lies in the detail that the sequence of (.) 
and |.| has been interchanged. The parameter $6 is nonzero if the tilt direction is locked 
to the nearest neighbor, next nearest neighbor, or to an intermediate direction. However, it 
would still be zero in a special case of locked state, where the tilt jumps between nearest and 
next nearest neighbors. In order to distinguish such a state from one where the tilt direction 
is really oblivious to the hexagonal lattice, we have also evaluated the related parameter $ 12 
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i n b 



6n . , , . 

J= i k=i 



(13) 



The parameter $6 and $12 are shown in Fig. [12| for fixed temperature T = 0.5e/fce as a 
function of pressure. At this temperature, the monolayer is tilted at all pressures shown. 
Fig. |7] demonstrates that the tilt direction is locked to the hexagonal lattice at low pressures, 
but apparently unlocks at II = 40e/<r 2 . That unlocked phases should exist in tilted hexatic 
liquid crystal films has been claimed by Selinger and Nelson ||. In crystalline phases, they 
are supposedly suppressed by the elastic interactions. Since our systems are too small to 
allow for a distinction between hexatic and crystalline order, they are obviously also too 
small to allow to decide whether the unlocked state is real or a finite-size artefact. 

In order to study the role of the chain flexibility, we have also performed a few shorter 
simulation runs (35.000 MCS) of systems with stiffer chains | 28fl . To this end, the stiffness 



constant k^ (cf. eqn. (g)) was increased by a factor often, k^ = lOOe. The area per molecule 
A/n, the melting order parameter \l/ 6 and the order parameter of collective tilt R xy for these 
systems are shown as a function of temperature for three different pressures II = 10, 30 and 
40e/<7 2 in Fig. 13. Up to the highest pressure II = 40e/<r 2 , the melting transition and the 



tilting transition are coupled. Moreover, the melting transition is shifted to much higher 
temperatures. This demonstrates once more that the melting transition in the system is 
basically driven by the chains. 



Our results for flexible chains are summarized in the phase diagrams Fig. O and Fig 



15| . We find at least four phases: the disordered fluid, an untilted ordered phase, two tilted 
ordered phases with tilt towards nearest neighbors and next nearest neighbors, and possibly 
an unlocked tilted phase. The areas per molecule of the two locked tilted phases are almost 
equal at the transition, even at low temperatures where the latter is strongly first order. At 
higher temperatures, the transition is so washed out that it cannot be located any more. 
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The transition between the tilted and the untilted ordered phase seems continuous. Between 
the tilted ordered phase and the disordered phase, it is presumably first order. The order 
of the transition between the untilted ordered phase and the disordered phase could not be 
determined, as discussed above. It should be stressed that none of our assertions on the 
order of the transitions has been corroborated by a finite size analysis, hence they should 
be regarded with caution. 

At surface areas per molecules smaller than A ps 0.8a 2 , i.e., at high pressures and low 
temperatures, the chains are squeezed together so closely that they form "rippled" structures 
where the beads of chains in neighbor rows are displaced with respect to each other in the 
z direction. This effect is clearly an artefact of our model and has not been investigated 
in detail, nor included in the phase diagram Fig. [13|. In the limit of vanishing pressure, 
on the other hand, the system has to assume a gas phase at all temperatures for entropic 
reasons. The transition between the gas phase and the condensed phase is subject to strong 
hysteresis effects at low temperatures. Nevertheless, we have been able to determine the area 
per molecule of the coexisting condensed state without too much computational effort on the 
basis of the following consideration: An upper limit is given by the area per molecule of the 
metastable condensed state at zero pressure, which does not decay within the simulation 
time at temperatures below T = 1.35e/ ks- A lower limit is provided by the area per 
molecule at the smallest pressure for which the transition temperature from the ordered 
to the disordered state has been determined, in our case II = le/a 2 . Since the areas per 
molecule do not depend strongly on the pressure in the condensed state, the coexistence line 
can thus be located fairly accurately (see Fig. |TJ|). 

Within the region of the disordered fluid, we have not found evidence for an additional 
liquid/gas transition. Such a transition would be expected at areas per molecule much 
larger than ~ 3cr 2 (where the critical point is found in two dimensional Lennard- Jones fluids 
f38fl), and correspondingly low surface pressures. We have spent some time searching for 
it, varying the temperature at very low pressure IT = 0.05e/cx 2 , and driving the pressure to 
zero at the temperature T = lAde/ks P5| . In a region around (II = 0.05e/<r 2 ,T w 1.7e//ce) 
or (II ps 0.04e/cr 2 ,T = lAbe/ks), the area per molecule varied rapidly, and strong density 
fluctuations were encountered. This suggests that liquid-gas critical point may be nearby. 
However, we have not been able to locate it so far. It may be hidden in the coexistence 
region. 



IV. CONCLUSIONS 

To summarize, we have studied in detail the phase behavior of a model of grafted 
Lennard- Jones chains, which is of interest as a "minimal" model for amphiphile monolayers. 
The model was found to show an impressive variety of phases, and its analysis gives useful 
insight into the mechanisms which drive some of the phase transitions in amphiphilic layers. 
In particular, it exhibits a disordered phase, an untilted ordered phase, and a number of 
tilted ordered phases, which are also found experimentally in Langmuir monolayers. The 
sequence of tilting transitions with increasing pressure (tilt towards nearest neighbors, tilt 
towards next nearest neighbors, no tilt) agrees with experiments and with earlier theoretical 
predictions. Furthermore, we have discussed the transition to the fluid state, and concluded 
from the form of the pair correlation functions in the different phases, and from the way the 
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transition temperature depends on the chain stiffness, that the transition is mainly driven by 
the chains, again in agreement with experimental [16] and theoretical [10,11,17] observations. 

In other respect, the phase diagram is still quite different from the experimental one (Fig. 
[I]). Some of the discrepancies are not surprising; for example, the model with its rotationally 
symmetric chains was never designed to reproduce the herringbone-ordered low-temperature 
structures. Other differences are more interesting. The pressure at the transition from 
the tilted to the untilted phase decreases strongly with temperature, whereas it is almost 
independent of the temperature in experimental systems. Likewise, the transition pressure of 
the swiveling transition between nearest neighbor tilt and next nearest neighbor tilt increases 
with temperature, whereas the line separating the Ov and L2 phase in Fig. [I] moves to lower 
pressures. This is presumably a consequence of the treatment of the head groups - more 
specifically, of the rigid constraints which are imposed on them in the model. The hard 
core interactions are much harder than the effective interactions between real head groups 
in water. Moreover, the heads in our model are confined to lie in a plane, whereas they can 
move in and out of the surface in real systems |40[| . 

Further refinements of the model will thus have to focus on the representation of the head 
groups. We have already mentioned the interplay between head size, spreading pressure, 
and tilting transitions. A more detailed study of the influence of the head size on the phase 
behavior shall be presented elsewhere [Q. Future work will be concerned with the effect 



of relaxing some of the constraints on the head groups, i.e., giving them additional degrees 
of freedom in the ^-direction, and possibly softening the interactions between them. One 
could also think of introducing interactions between the tails and the substrate. However, 
the tails hardly come into contact with the substrate at most densities of interest, therefore 
this will probably not change the phase behavior significantly. 

On the other hand, we have seen that already the present simple model reproduces many 
important properties of amphiphile monolayers. Hence it can be used as a starting point for 
further investigations. In particular, simulations of much larger systems and a systematic 
variation of system sizes would be desirable to shed light on some of the questions which 
have remained open in the present study. These would help to elucidate the exact nature 
of the tilting transitions and the order of the melting transition, to examine the unlocked 
tilted state, and to clarify whether our model actually does exhibit hexatic phases. 
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Temperature 



FIG. 1. 

Generic phase diagram for fatty acid monolayers. The phases LS, S and CS are on average 
untilted, whereas Ov and L' 2 show tilt towards next nearest neighbors, and L 2 , L 2 towards 
nearest neighbors. In CS, S, L 2 and L 2 , the backbones of the hydrocarbon chains are 
ordered. In CS and L 2 , the molecules have crystalline order in addition. 
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FIG. 2. 

Area per molecule A/n in units of a 2 vs. temperature T in units of e/ks for a choice of low 
(a) and high (b) pressures II (in units of e/cr 2 ) as indicated. 
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FIG. 3. 

Radial pair correlation functions g(r) vs. r in units of a at pressure II = lOOe/cr 2 and 
temperature T = le/k B - Correlation functions are shown for the heads (solid line), for the 
points where the molecules cross the plane at z = 2a above the surface (dotted line), and 
for the projection of the center of gravity onto the xy plane (dashed line). The values of 
g{r) for T = Q.le/k B are divided by a factor of five for the clarity of presentation. 



17 



(a) 



4 
3.5 

3 
2.5 

2 
1.5 

1 

0.5 




1/5 *T = 0.10 
T = 1 .40 
T = 1.41 
T = 2.00 




(b) 



9 



(r) 



4 
3.5 

3 
2.5 

2 
1.5 

1 

0.5 




1/5 *T = 0.10 
T = 1 .40 
T = 1.41 
T = 2.00 




FIG. 4. 

Radial pair correlation functions g(r) vs. r in units of a at pressure II = le/a 2 and various 
temperatures as indicated. Correlation functions are shown for the heads (a), and for pro- 
jections into the xy plane of the centers of gravity (b). Temperatures T are given in units 
of e. The correlation functions g(r) for the temperature T = O.le have been divided by a 
factor of 5. 
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FIG. 5. 

Structure factor S(q) in the xy plane (q z = 0) for a disordered state (a) and an untilted 
ordered state (b). Parameters are II = 10e/cr 2 , T = 2.5e/k B in (a), and II = 50e/a 2 , 
T = 2.0e/k B in (b). 
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FIG. 7. 

Order parameter \l/ 6 vs. temperature T in units of e/ks for different pressures II (in units 
of e/a 2 ) as indicated. 
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FIG. 8. 

Order parameter R xy in units of a 2 vs. temperature T in units of e/ks for different pressures 
II (in units of e/<r 2 ) as indicated. 
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FIG. 9. 

Average tilt angle (9) in degrees vs. temperature T in units of e/ks for different pressures 
II (in units of e/o" 2 ) as indicated. 
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FIG. 10. 

Product of the cosine of the average tilt angle with the area per molecule, Acos({9))/n, in 
units of a 2 , vs. temperature T in units of e/ks for different pressures II (in units of e/o" 2 ) 
as indicated. The horizontal line indicates the position of a c = 0.985(j 2 . 
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FIG. 11. 

Histogram P{6)/ sin# of the tilt angle 6 (in degrees) at pressure II = 50e/cr 2 for different 
temperatures (in units of e/fce). 
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FIG. 12. 

Order parameters $6 an d $12 vs. pressure II in units of e/a 2 at temperature T = 0.5e//ce. 
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FIG. 13. 

Order parameters and R xy , and area per molecule A/n vs. temperature T in units 
of e/ks for pressures II = 10e/a 2 (filled circles), 30e/a 2 (open squares), 40e/cr 2 (stars) in 
systems of stiff chains, (kA = lOOe). 
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FIG. 14. 

Phase diagram in the pressure-temperature plane. Pressure n is given in units of e/cr 2 , 
and temperature T in units of e/k B . LE denotes disordered phase, LC-NN ordered phase 
with tilt towards nearest neighbors, LC-NNN ordered phase with tilt towards next nearest 
neighbors, and LC-U untilted ordered phase. The transition between LC-NN and LC-NNN 
could not be located at pressures above II = 20e/a 2 . See text for more explanation. 



28 




Area per molecule A/n 



FIG. 15. 

Phase diagram in the area-temperature plane. Area per molecule A/n is given in units 
of a 2 , and temperature T in units of e/ks- LE denotes disordered phase, LC-NN ordered 
phase with tilt towards nearest neighbors, LC-NNN ordered phase with tilt towards next 
nearest neighbors, and LC-U untilted ordered phase. See text for more explanation. 



29 



